Electron-vibration interaction in transport through atomic gold wires 
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We calculate the effect of electron-vibration coupling on conduction through atomic gold wires, 
which was measured in the experiments of Agrait et al. [Phys. Rev. Lett. 88, 216803 (2002)]. The 
vibrational modes, the coupling constants, and the inelastic transport are all calculated using a 
tight-binding parametrization and the non-equilibrium Green function formalism. The electron- 
vibration coupling gives rise to small drops in the conductance at voltages corresponding to energies 
of some of the vibrational modes. We study systematically how the position and height of these 
steps vary as a linear wire is stretched and more atoms are added to it, and find a good agreement 
with the experiments. We also consider two different types of geometries, which are found to yield 
qualitatively similar results. In contrast to previous calculations, we find that typically there are 
several close-lying drops due to different longitudinal modes. In the experiments, only a single 
drop is usually visible, but its width is too large to be accounted for by temperature. Therefore, 
to explain the experimental results, we find it necessary to introduce a finite broadening to the 
vibrational modes, which makes the separate drops merge into a single, wide one. In addition, we 
predict how the signatures of vibrational modes in the conductance curves differ between linear and 
zigzag-type wires. 

PACS numbers: 72.10.Di, 73.23.-b, 73.40.Jn 
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I. INTRODUCTION 

The electron transport properties of atomic point con- 
tacts between two metallic electrodes have been inten- 
sively studied during the past decade 0. Contacts of 
this type are typically formed by using mechanically con- 
trollable break junctions (MCBJ) or with the tip of a 
scanning tunneling microscope (STM) . It has been found 
that the conductance of such contacts depends strongly 
on the electronic structure of the metals, and for monova- 
lent metals there is a tendency for quantization in units of 
the conductance quantum Go = 2e^//i 2]. Point contacts 
formed from Gold (Au), Platinum (Pt), or Iridium (Ir) 
by one of the MCBJ or STM methods have the further 
interesting property of being able to sustain single-atom 
thick chains, so called atomic wires Following 
their discoveryja good amount of experimental 0,0 and 
theoretical U 0, [3> ^> ^> El work has been carried 
out to further investigate the conduction properties of 
atomic wires. It is by now well established that the zero- 
bias conductance of gold wires is close to one Go due to 
a single, almost fully o pen transmission channel at the 
Fermi energy [2, HI, llM Ea ■ However, less detailed work 
has been done in the study of truly nonequilibrium prop- 
erties, such as the current-voltage characteristics 17] . 

In recent experiments, the conductance vs. voltage 
characteristics G{V) of gold wires formed by the STM 
technique were measured 18]. It was observed that the 
conductance often has a very pronounced, single drop 
from Go at a critical voltage Vph = 10 — 20 mV, mark- 
ing the onset of a dissipative process. The size of the 
drop was on the order of 0.5% — 2.0% of Gq. It was 
also found that stretching of the wire typically leads to 



an increase in the step, and to a decrease in the criti- 
cal voltage Vph- Based on simple arguments for infinite 
single-orbital tight-binding chains, these findings were in- 
terpreted as a sign of the excitation of vibrational modes 
in the wire: only a single longitudinal mode with twice 
the Fermi wave vector can be excited, since this corre- 
sponds to the momentum which must be transferred from 
an electron to the vibrations in a single backscattering 
process. Although the validity of such arguments for a 
wire of finite length (of typically less than 10 atoms) can 
be questioned, the interpretation was backed up by first- 
principles calculations |1^ l20l ] . The authors of Ref. 
emphasize the importance of so-called alternating bond 
length (ABL) modes, and in particular the longitudinal 
mode of highest frequency. 

Although it seems evident that the interpretation 
based on vibrational modes is essentially correct, what 
is still lacking is a systematic study of the behavior of 
wires with varying numbers of atoms, and surrounded 
by various lead geometries. Many questions of the basic 
physics are also still not very well understood: When ex- 
actly does the electron-vibration coupling lead to a drop 
and when an increase in the conductance? Why does 
there appear to be just a single drop in the experiments 
of Ref. list although the momentum conservation is not 
exact? What determines the height and width of this 
drop? Below we aim to discuss the possible answers to 
some of these questions. 

In this paper we concentrate on studying the current- 
voltage characteristics of gold wires. We use a Slater- 
Koster JlJ type tight-binding (TB) approach, where 
the parameters are taken from the non-orthogonal 
parametrization of Papaconstantopoulos and coworkers 
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11^ 11^ l^sf . The use of such a parametrization 
[2a. 1271 l2q makes the modeUng of atomic wires com- 
putationally less intensive as compared with fully ab ini- 
tio methods. The approach is still microscopic in that 
it takes into account the symmetries of the atomic s, p, 
and d valence orbitals, which, via hybridization, form the 
conduction channels. It is also general enough to allow 
one to model everything within the same framework: we 
use the parameters to compute the total energy of the 
wire, and thus optimize the geometry. After this, the nor- 
mal modes of oscillation and the electron-vibration cou- 
pling constants can be computed. Finally, we calculate 
the transport properties using the non-equilibrium Green 
function (or Keldysh) approach. Our implementation is 
very similar to that of Ref. 28, and the present work is, 
in essence, an extension of that to inelastic transport. In 
addition to the full ab initio calculations 19], the effect 
of electron-vibration interactions on transport through 
molecular wires has been rec ently studied by some sim- 
ple single-level models |2^|33,|31|- The tight-binding ap- 
proach stands somewhere in between these two extremes. 

We compute the conductance to lowest nontrivial order 
in the electron- vibration coupling constant. There are es- 
sentially two well-defined limits which can be studied. In 
the first limit the vibrational-mode distribution remains 
in equilibrium due to a strong coupling of the modes to 
an external equilibrium bath formed by the leads. In the 
opposite limit the distribution is driven to strong non- 
equilibrium by the bias voltage. These are the extern ally 
damped and the externally undamped limits of Ref. Il9l 
However, in the first limit one should also account for the 
strong broadening of the vibrational modes. We derive 
equations which take this into account in a phenomeno- 
logical manner. 

With our simple, self-contained method of optimizing 
the geometry, we obtain vibrational frequencies which 
are of the correct order of magnitude, usually to within 
a factor of two. We study how the positions and heights 
of the conductance drops due to the electron-vibration 
coupling vary with elongation of a linear wire, and find 
a g ood overall agreement with the experiments of Ref. 
Il8l The height of the conductance steps grows together 
with the length of the wire, being typically of the or- 
der 0.5% — 5% of Go for wires of 11 atoms or less. 
As found in the earlier calculations we find that 

the highest-frequency longitudinal modes usually couple 
most strongly, although there seems to be no fundamen- 
tal reason for a bias toward the "ABL" modes. But in 
contrast to previous theoretical results, the conductance 
drop is usually found to occur in two or more consecutive 
steps which are due to several close-lying longitudinal vi- 
brational modes. Thus we find the "mode selectivity" 
to be only very approximate. However, the steps can 
be made to merge into a single one, when we introduce 
a large enough phenomenological broadening to the vi- 
brational modes, such that the experimentally observed 
step widths of ~ 5 meV are accounted for. We also briefly 
study the conductance signatures of chains which have a 




FIG. 1: Geometry A, without "pyramids". A zigzag wire with 
Nch ~ 6 atoms is shown. 




FIG. 2: Geometry B, with "pyramids" and a linear wire of 
Nch ~ 3 atoms. 

zigzag-like form, instead of the linear one. 

The paper is divided into the following parts. In Sec. 
im we start by defining the problem, and discussing the 
electron- vibration coupling. In Sec. lIIII we briefly discuss 
the calculation of the vibrational modes and the electron- 
vibration coupling constants, as well as our methods of 
computing the transport. After this. Sec. IIVI introduces 
the important wide-band approximation to the full for- 
malism. In Sec.^we use the formalism to analyze simple 
tight-binding models, and in Sec. IVII the full spd-tight- 
binding model. Section IVIII ends with some conclusions 
and discussion. Technical details which are not of im- 
mediate importance are postponed to the appendixes. 
These include a discussion of the calculation of the ma- 
trix elements needed for the coupling constants. Readers 
mainly interested in the results can directly jump to Sees. 

Ellvnl 

II. DEFINITION OF THE PROBLEM 

To model transport through atomic wires, we consider 
two idealized geometries, shown in Figs.^andl^l We call 
these geometry A and geometry B, respectively. Both 
involve a gold chain of Nch atoms suspended between 
two gold leads. As the leads, we simply use semi-infinite 
"bars" , where the repeat unit consists of two layers, with 
12 and 13 atoms respectively, mimicking an infinite fee 
[001] surface, where the z axis is always chosen parallel to 
the axis of the wire. The particular choice for the leads 
should not be very important, as long as they are infinite 
in one direction, and wider than the contact region. In 
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geometry B, the chain connects to small clusters of atoms 
or "pyramids" on the surfaces (in our case consisting of 
9 atoms), making it perhaps the more realistic one of 
the two. For technical reasons, the geometry is divided 
into three parts, the semi-infinite left (L) and right (R) 
leads, and the "central cluster" (C), which also includes 
the pyramids if any. These parts are indicated in the 
figures. We shall also consider some simple models where 
the regions L, C, and R are parts of an infinite linear 
wire. 

Our objective is to model the effect of vibrations (or 
"phonons") of the wire on the transport, when a volt- 
age is applied over the contact. Within a tight-binding 
picture, the system of electrons coupled to vibrational 
modes is described by the Hamiltonian 

H = He + Hyib + He-vib, 

where 

ij 

Hvib = ^ fl'^abl.ba (^ l^ 

a 

ij a 

Here Wq, are the vibrational frequencies, Hij — {i\H\j) 
are the matrix elements of the equilibrium single-electron 
Hamiltonian H in the atomic-orbital basis and Xfj 

are the electron-vibration coupling constants. The in- 
dex i denotes collectively the atomic sites and orbitals, 
and a runs from f to 3Nyii,, where Nyu, is the number 
of atoms in the system, which are allowed to vibrate. 
The creation and annihilation operators for vibrational 
modes foJ,,6Q satisfy the bosonic commutation relation 
[foajfe^j] = 6a/3. The electronic basis is in general non- 
orthogonal, with overlap matrix elements Sij = {i\j). 
Thus the anticommutator for electron operators is 
given by = [S'%. 

Hereafter we denote the matrices with components Aij 
with a boldface symbol A. The matrices H, S, and A" 
are all symmetric in our case. In the spd TB model, 
the matrix elements Hij and Sij are obtained directly 
from the parametrization l^^. These can also be used 
to calculate the the vibrational frequencies uja and the 
coupling constants A", as we shall now describe. 

III. METHODS 

The solution of the inelastic transport problem involves 
a few rather separated sub-problems: the optimization 
of the geometry and evaluation of the vibrational modes, 
estimation of the electron- vibration coupling constants, 
and finally the calculation of the transport. In the fol- 
lowing we give only a brief description of each of them, 
and refer the reader to appendixes for details. Our basic 



approach is to solve for the elastic transmission prob- 
lem exactly, and then to take the electron-vibration cou- 
pling into account in a slightly modified version of lowest- 
order perturbation theory. Other works have considered 
the so-called self-consistent Born approximation (SCBA) 
|l9l l3l| , where some of the terms in the perturbation ex- 
pansions are effectively summed to infinite order. How- 
ever, this is not essential for describing the basic physics 
which is involved in the present problem. 

A. Vibrational modes and the electron-vibration 
coupling constants 

The calculation of the vibrational modes requires 
knowledge of the total (ground-state) energy of the sys- 
tem as a function E{Rk) of the ionic coordinates Rk with 
A: = 1, . . . , Nyii,. This energy must be minimized to find 
the equilibrium configuration R^^'^ . Now consider small 

displacements Qfc = Rk — R^*^ around the equilibrium. 
The Hamiltonian describing the oscillations of the ions 
around i?^"' is given by 

where Mk are the are the ionic masses, ^,1^ — x,y, z de- 
note Cartesian components of vectors and H is the Hes- 
sian matrix: Hkfj.,iL> = d^E/dRkfj.dRi^. This can be diag- 
onalized by the transformation — J2a=i'' ^kfi.ala, 
where Qa are the normal coordinates. Thus, we obtain 

Hion = \ Y^aiia + '^a^a); whcrC LUa (a = 1, . . . , 37V„i;,) 

are the vibrational frequencies. The transformation ma- 
trix A is normalized according to A^ MA = 1, M being 
the mass matrix — in our case M is simply a scalar 
giving the mass of a gold atom. Upon using the canoni- 
cal quantization prescription qa — {hjlua^^"^ (6^, + 60), 

qa ~ i {huJa/'2,y^^ {bl^ ~ ba), one finally obtains Hvib m 
Eq. O. 

The electron-vibration interaction may be derived as 
follows Is^i ' Assume that the electronic single- 
particle Hamiltonian 77 is a function of the ionic coor- 
dinates, denoted collectively as R. Then we may expand 
i?(i?(o) + g) « iJ(i?(o)) + EfeQfc • Vfei?|^^o. Defining 

H'e = J2ij d-l{i\H{R^'^'> +Q)\j)dj , inserting the expansion, 
and using the canonical quantization for qa again, one 
finds H'e = He+He-vrb [cf. Eq. Q], where H = i/(i?(°)), 
and the electron-vibration coupling constants are given 

by 




where M^J" — («| Vfe^i/|jj^g|j). The calculation of these 
matrix elements is explained in Appendix^ In Eq. ^ 
we have added a dimensionless factor Aq to describe the 
strength of the coupling — in the physical case Aq = 1. 



4 



B. Propagator formalism 

Use of a local basis allows one to partition the elec- 
tronic Hamiltonian and overlap matrices into parts ac- 
cording to the division in L, C and R regions: 







Hlc 






Sll 


Slc 


Slr 


H = 


HcL 


Hcc 


HcR 


, s = 


ScL 


Sec 


ScR 






Hrc 


Hrr 




Srl 


Src 


Srr 



Although the dimension of the problem is infinite, its 
single-particle nature allows for very effective methods of 
solution, as long as we may assume that H j^l — HJ^j^ « 
and Srl — S^r « 0, which we shall do. We shall use 
the method of non-equilibrium Green functions (NEGF) . 
In this method, one can restrict the problem only to the C 
part by introducing energy-dependent lead self-energies 
which take into account the presence of the semi-infinite 
L and R leads in an exact way. 

The quantity from which all elastic transport proper- 
ties may be extracted, is the retarded Green function 
of the C part in the absence of electron-vibration cou- 
pling. We call it G , and it may be written as G (e) = 
[eScc - Hcc - S£ - T,'r]-\ The lead self-energy S£ 
is given by = tcLQlLtLC, and Tl = i(S£ - S^), 
where we define tcL — Hcl ~ ^Scl- The matrix 
91l{^) = [(e + hL/2)SLL - Hll]-^ is the lead (sur- 
face) Green function, where 7^ = O'^. Similar equations 
hold for S^. The lead Green functions g^^^ and are 
"surface" Green functions for the semi-infinite leads. We 
compute these with the so-called decimation technique 
[3^ , using TB parameters for bulk in the case of the 
full spd model. The electron-vibration interaction gives 
rise to further self-energies, as will be discussed below. 

The vibrational modes should in principle be treated 
in an analogous way, by introducing lead self-energies for 
their propagators. However, here we restrict the modes 
strictly to the wire of Nch atoms within the C region (i.e. 
Nvib = Nch) and use the corresponding normal-mode 
basis for them. Thus the number of modes which we 
have to consider is only 3Nch, and their "lead coupling" 
is taken into account only in a phenomenological way. 

More details on the propagator technique, including 
the expressions for the phonon propagators and all self- 
energy diagrams, are given in Appendix lUl 



C. Calculation of current 

The most important physical observable which we are 
interested in is the electric (charge) current through the 
atomic wire, when a voltage V is applied. We denote 
eV = fJ-L^^ l^-R, where fiL.R. are the L and R side chemical 
potentials, and e > is the absolute value of electron 
charge. We also define /L,fl(e) = /(e — fJ-L.R): where 
/(e) = l/[exp(/3e) -I- 1] is the Fermi function, and /3 — 
l/ksT is the inverse temperature. 



It may be shown that the current flowing through the 
interface from LtoC{C to R) in stationary state is given 
by (see Appendix IB|| 

/"^-iy/ ^Tr[Gin{^)tnc{e)-tcn{e)G+ci^)], 

.(3) 

where 51 = L{R) is chosen with the upper (lower) sign, 
and the factor 2 accounts for spin degeneracy. The Green 
functions G^ are defined in Appendix lUl bv Eq. ljGl|l . 
Developing Eq. Q further, it is convenient to split it into 
two parts: /^'^ = lei + Xj„e; , where 

X [(/l,_R - 1)S+J„-^ - fL,R'^e-vib\}- 

(4) 

Here we define the full retarded and advanced Green 
functions G'^''', where C = [eScc ' Hcc -^l - ^r - 
K-vtb]~^ and G" = [GH^ The new self-energies K-v^b 
and S^J'^j^ are due to the electron-vibration interaction, 
and they are discussed in more detail in Appendix [CI 
Since they vanish in the absence of A", we call the I^^^i 
jjar t an "inelastic" current, while Igi is the "elastic" part 

If we do lowest-order perturbation theory with respect 
to A", we may expand G^ — G + G Sg_^j^G -I- • • • . 
In this way the elastic current is split into two parts as 
/pj = I^i + Slei , where Slei is an "elastic correction" . We 
find 

lei-^J ^Tr[G^rHG"r^](/z.-/i.) 
Slei = y 1 ^ReTr[riG''s:_„,,G''rflG"] 

X (/l - /r) 

X [(/l,_R - l)^^~vib ~ fL,R'^~Xib]}- 

(5) 

A proof of current conservation, that is = = I, is 
sketched in Appendix IdI 

Besides the charge current, other interesting observ- 
ables would be the heat current (or power dissipation) 
|l9l|. current noise js^ . and possibly spin current in case 
of magnetic materials. We only consider the charge cur- 
rent here, as it is the only one which can easily be mea- 
sured. More specifically, we shall be interested in the 
differential conductance G{V) = dl/dV and its deriva- 
tive, since these quantities reveal the signatures of the 
vibrational-mode coupling most clearly. 
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IV. WIDE-BAND LIMIT 

Even in the case of the perturbative current formu- 
las [Eqs. (0], the expressions will involve double en- 
ergy integrals which can be very cumbersome to evaluate. 
These general formulas are discussed more in Appendix 
IeI where they are rewritten in terms of distribution func- 
tions and energy-dependent transport coefficients. How- 
ever, the existence of different energy scales in the prob- 
lem allows us to make an important simplification. 

The energies of the vibrational modes are on the order 
of 10 meV, so that we are only interested in the dif- 
ferential conductance for voltages up to « 40 mV, 
at most. Together with the temperature T f« 4.2 K, 
this determines the width of the energy window around 
the Fermi energy {ep) which is important for transport. 
However, for the atomic wires which wc are considering, 
the electronic density of states tends to vary at the much 
larger energy scales ~ 1 eV around ep- Thus, to a good 
approximation, we may neglect this energy dependence, 
and simply evaluate all the electronic Green functions at 
the Fermi energy. This approximation is often called the 
"wide-band hmit" (WBL). 



Eq. linHl, that is 

( \ ^ 1 _ 1 ^/2 

P"^^' TT (e - fi£^„)2 + ry2/4 ^(e + fito„)2+,72/4' 

(8) 

where we take 77 as a finite phenomenological parameter 
describing the effect of coupling the vibrational modes to 
an external bath. This bath is provided by the leads [sif . 
However, we are neglecting any renormalizations of the 
bare frequencies oja, so that the main purpose of 77 here 
is to broaden the DOS. Similarly, Redjj is obtained from 
Eq. iHSll . 

The current involves a number of transport coefficients, 
which determine the shape of the current-voltage char- 
acteristics. The coefficients of py™ may be computed as 
follows: 

T,, = Tv[CfTRG''TL]\,,, (9) 



T^" = Tr[G'rj^G"A"G°riG'A"] , (10) 



= 2RcTr[G''rflG°rLG'^A"G''A"]|,,, (11) 



A. Current 

In the WBL approximation, the current expressions of 
Appendix ^ may simplified considerably, since some of 
the energy integrals may be done analytically. In this 
case the current / is easily divided into symmetric and 
asymmetric parts, according to the symmetry of their 
contributions to G{y) under the reversal of V 3^. Thus 
J _ jsyrn ^ jasy ^ .y^f}jg]-g ^}jg symmetric Current is 



rj 2 o POO 

jsym ^ _j.^y + V / dc.ip„(^i) 

[(T- + T™)(2iV„(c^i) + 1) - (T-^« + 7:")]ey 

ecLR . rr.in^ ( UJl - eV Ui + eV ' 



-(t; 



(6) 



and the asymmetric part is 
2e 



r 



asy _ 



h 

2n(uJi)LjJi 



/OO 1 



uii — eV 



(jji + eV 



(7) 



Here n(e) = l/[exp(/3e) — 1] is the Bose distribution and 
Na is the voltage-dependent mode distribution, to be dis- 
cussed shortly. The function is the vibrational density 
of states (DOS), given in general by Eq. (|C12|) . We ap- 
proximate it here by using the imaginary part of 0?^ in 



rpecLR rpecL i rpecR 

^ReTT[G"TRG"TLG'\"iG' - G")A"] 



(12) 



Here T^'^^^ is a sum of the two coefficients T^'^^ and 
T^'^"- defined in Eqs. The coefficient T™ is always 

positive, while T^^^^ and T^'^ can apparently have ei- 
ther sign. The latter two represent interferences between 
various processes, and are responsible for the enhanced 
backscattering needed for the conductance drops. The 
asymmetric coefficient 



T^^y = ReTr[G'r;?G°rLG''A"G''(r;^ - Tl)G''\'']\,, 

(13) 

vanishes for sufficiently symmetric junctions, making 
jasy _ Q us, this is always the case, even with 

zigzag chains. 



B. Mode distribution 

Our approximation is essentially that of lowest-order 
perturbation theory in the coupling constant A". How- 
ever, as explained in Appendix there is a natural 
way of extending the theory somewhat further by tak- 
ing into account the "phonon polarizations" in a voltage- 
dependent distribution function Na ■ This is given by Eq. 
(ITTll . i.e. 



l lmn+-{t) + n{e)r]e/nuja 
'2 ImHj;(e) - r;e/2fta;a 



(14) 
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where 77 is the same bath-coupHng parameter as in pa- 
Here Imn+~ and Imll^ are imaginary parts of the 
phonon polarizations (see Appendix Q • In the WBL 
one may show that 



27rlmn+"(e) 

e + ey e-eV 



(15) 



and 



Imn;(e) « -(e/7r)Tr[A"ImG A"ImG ], (16) 

which is proportional to the electron-hole damping rate 
of Ref.EE Using G\TL+TR)Cf = -2ImG\ one may 
easily show that ior V — Eqs. 1)14(1 - ^6(1 indeed yield 
Na{e) — n{e) for any 77. 

Here it is important to note a few things. In the ex- 
pression for the distribution function, the limit 77 — » 0-1- 
corresponds to the case where the vibrational modes are 
uncoupled from leads. Supposing that one also wishes 
to take the phonon polarizations to zero, which is for- 
mally accomplished by taking Aq ^ 0, one discovers an 
interesting thing: the two limits do not commute. 

If we take first the limit 77 ^ 0, then the result ac- 
tually becomes independent of Aq, since HJj,H+~ cx Aq 
and the Ag-factors cancel. A physical interpretation can 
be described as follows. If the vibrational modes are not 
coupled to any external bath, then even an infinitesimally 
small coupling constant can eventually lead to a station- 
ary state with a strongly nonequilibrium mode occupa- 
tion. Here emission and absorption of phonons are in 
balance, and hence there is no net energy transfer be- 
tween the electrons and the vibrations. Following Ref. 
[T9I we call this the externally undamped limit, although 
our way of computing Na is quite different. In this case 
the voltage-dependence of NaifkOa) shows a sharp kink 
at V = huja/e, and a subsequent linear increase [13 • 

In the opposite case, where Aq — > first, the expres- 
sion becomes independent of rj, and we recover the Bose 
distribution. This corresponds to the limit where the vi- 
brational modes are strongly damped by coupling to a 
heat bath which is in equilibrium. This is the externally 
damped limit. However, for a finite A" this limit can 
only be reached with a large enough finite rj. Thus the 
externally damped limit should also imply a considerable 
broadening of the vibrational modes. 

Note that in both of the above limits, Na is zeroth 
order in Ao. In these two cases our expression for I'^v"^ 
[Eq. jnj] is indeed of second order in Aq , and there should 
be no corrections within that order. In general, however, 
Eq. I(14() for Na generates terms of all orders in Aq. These 



are not, strictly speaking, warranted, because they repre- 
sent only a small subset of all possible higher-order terms 
in the current. 



C. Further discussion 

Apart from the pct-weighted integral, the only differ- 
ence between Eq. and the approximation of Ref. 
for ry"" is that T;^" T^<^lr^ rpj^g difference T^^-T^"^^ 

~ r 

is proportional to Re G , and is typically very small. It 
could well be neglected. If Na =constant, then the di- 
rection and size of the conductance step is solely deter- 
mined by the combination t;^<'LR^ j.vn_ q^^^^ j.m 

positive, the Ii„ei part of the current always tends to 
increase the conductance. This can be seen as a result 
of the appearance a new, inelastic conduction "channel" 
when eV > hiOa- However, numerically one finds that the 
coefficient Ta'^^^, due to interferences between various 
elastic processes, is typically negative and \Ta''^^\ > 
when To R:! 1. Thus, the net effect of the vibrational 
coupling is a decrease in the conductance. In general, 
the voltage-dependence of Na also affects the shape of 
G{V). Since T^" + is usually also negative, the ef- 
fect of local heating is to give a finite negative slope to 
G{V) after the step. There can also be an increase in the 
apparent height of the step. 

By using the vibrational DOS of Eq. ||HJ), we are ne- 
glecting broadenings and frequency shifts due to the 
electron-vibration coupling. The shifts are given by the 
quantities ReH^ and, using Eq. IjClOp . we estimate them 
to be on the order of —1 meV. Neglecting the broadenings 
is certainly justified, since they are typically on the order 
of \lmlla{fiuja) \ / fioJa ^ 10~^, which is much smaller than 
the usual kBT/fkUa due to temperatures T ~ 1 K |19| . 
If we assume also 77 to be small compared to the other 
relevant energy scales, i.e. if |ImH^|, 7; <C fc^T ^ ftw^, 
then Pq(lji) « S{lui — fkUa) — S{uJi + huja), and the uji 
integral in Eq. lO may be done analytically. The calcula- 
tion of the lead-coupling effects is difficult, but it is very 
possible that they can lead to broadenings rj/fkUa 1. 
Thus, the validity of the delta-function approximation of 
Pa may be questioned. Naturally, the lead coupling can 
also shift the vibrational frequencies. 

Note also that Na{e) diverges if /3e 0. Since in 
Eq. lO Na{e) is evaluated at e ~ huJa, there can be an 
arbitrarily large renormalization of the zero-bias conduc- 
tance if there are very low-frequency vibrational modes 
{huJa ksT) with a strong coupling to the electrons. In 
such a case, the theory appears to break down. With 
linear wires the condition ksT <^ fkUa is easily satis- 
fied for all strongly coupled modes a, but in case of 
the zigzag-chains (and/or other materials) it may not 
be. In any case, we limit our study mostly to the linear 
wires in this paper. We note that these restrictions are 
present also in Ref. [s^, where 7; = 0"^. In fact, in our 
formulation the problem is perhaps partly corrected by 
the presence of the "broadening factor" pa- In general 
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FIG. 3: Infinite multi-orbital nearest-neighbor chain in equi- 
librium, where the sites are separated by a distance a. 



dujipa{(^i) < 1, and when huja — > 0, the whole func- 
tion tends to zero. Thus low-frequency modes contribute 
to the current with a very small weight. 



V. SIMPLE CHAIN MODELS 

Here we shall first present some example results, rising 
an adaptation of the simple chain models of Refs. 1371 and 
l20l Only after this we turn to the full spd-tight bind- 
ing parametrization. Chain models of this kind are very 
appealing, because they can be studied analytically to a 
large extent, and allow us to make our main points in a 
simple fashion. 

The single-particle Hamiltonians H which we consider, 
are of the block-tridiagonal form 



H 



(17) 



graphically depicted in Fig. |3| Here e includes the on- 
site energies, and are the inter-site hopping matri- 
ces between sites i and i -I- 1. They are modulated as a 
function of the longitudinal atomic displacements Qi as 
ii,i+i{Q) = i° +t'{Qi — Qi+i)- We only consider orthog- 
onal tight-binding here, such that S = 1. The chain is 
split into three parts, where the C part has Nch atoms. 
The L and R "leads" are semi-infinite chains. The dis- 
placements Qi are restricted to the C part, where the 
Hessian is assumed to be of the form 



-1 2 



(18) 



with fixed boundary conditions at the ends [23 • Unlike 
in the spd model to be explained below, here the spring 
constant K is taken as a new, separate parameter. The 
vibrational modes are simply obtained by diagonalizing 

this n. 



A. Single (s) orbital: discussion 

To remind ourselves of some of the well-known argu- 
ments, let us first discuss the simple single-orbital or- 
thogonal model of Ref . 0| . 



If we assume that charge neutrality is achieved by one 
electron per atomic site, then the Fermi energy ep lies 
exactly in the middle of the band — 2 1 | < e — eo < 2 1 to | , 
where eq is the on-site energy, and to is the inter-site 
hopping. The Fermi wavevector is then kp — 7r/2a, 
where a is the inter-atomic distance. Since the chain 
has full translational invariance, momentum conservation 
dictates that only vibrational modes with the wave vector 
Qyib = 2kp — t: / a may be excited, via the backscattering 
of Fermi-point electrons. Furthermore, if the atomic or- 
bitals are invariant with respect to rotations around the 
axis of the chain (s-wave, say) the modes can only be 
longitudinal ones. As it happens, qmb =1^/0, corresponds 
exactly to the Brillouin zone boundary and thus to the 
highest-frequency modes. 

More technically, the momentum conservation can be 
seen to follow from the form of the Fourier-transforms 
Mi^ fe2 ^ki-k2,q of the coupling matrix elements M^j 
|32| . We note that the momentum conservation can only 
be strictly valid if one considers the actual vibrational 
modes of the whole infinite chain. Interestingly, it ap- 
pears to remain approximately valid even if the vibra- 
tional modes are restricted only to a small, finite part 
of the chain Indeed, as found in Ref. I20I when the 
Fermi energy is in the middle of the band, there is only 
a single visible step in the conductance. But, as we shall 
discuss next, this appears to be somewhat accidental. 



B. Two (s and pz) orbitals 

Although such a single-orbital model can correctly de- 
scribe the most important experimental observations, it 
neglects some details of realistic gold wires. Although 
atomic gold contacts and chains typically have only one 
fully open conduction channel at the Fermi surface, this 
channel is actually formed from the hybridization on mul- 
tiple orbitals with rotational symmetry around the chain 
axis (s, Pz, d^z^-r'^). Thus we consider here a simple 
generalization of the above chain model to the case of 
two orbitals, with s and Pz characters, respectively. In 
Table ^ we show a set of example parameters for this 
model. The resulting transmission and DOS curves are 
shown in Fig. 0] and the conductance-voltage charac- 
teristics in Fig. 13 The density of states is defined as 
I?i(e) — — (I/tt) ImG'M_iQ(e), where i now stands for 
atomic sites and a for orbitals. 

The results are very similar to those of the single- 
orbital chain '2^, but the situation is a bit closer to 
what happens in the full spd model. For example, the 
hybridization of s and pz result in a gap in the DOS. If 
we still populate the chain with one electron per site, the 
Fermi energy will again be in the middle of the lower band 
(the "s band"), and there is only a single drop in the con- 
ductance. This is illustrated by the ep = —0.4 eV case 
in Fig. 1^1 But for ep's deviating from the center of the 
band, we find that the conductance step generally con- 
sists of several substeps, although the total step height 
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Quantity 



Symbol Value 



^SS 



+0 







•'sp 



Number of chain atoms 
s-orbital energy 
p-orbital energy 
Bare s — s hopping 
Bare p — p hopping 
Bare s — p hopping 
s — s hopping modulation t'^^ 
p — p hopping modulation t'^^ 
s — p hopping modulation t'^p — 
Fermi energy ep 
Atomic Mass M 
Spring constant K 
Temperature T 
Phonon broadening rj 



0.0 eV 
1.0 eV 
-0.5 eV 
0.3 eV 
0.35 eV 
-0.3 eV 
0.3 eV 
0.3 eV 
-0.4—0.2 eV 
197 a.m.u. 
2.0 eV/A^ 
1.0—4.2 K 
0.002 meV 



TABLE I: Parameters for the sp^ schain. Here t^j are elements 
of the matrix iP, for example. 




FIG. 4: Elastic transmission To (solid line) and density of 
states (dashed line) for an spz chain model. 



remains almost the same. Notice that if we were to use 
two electrons per site, then the Fermi energy would lie 
in the gap, and the transmission would be zero. It must 
be noted, however, that the assumption of the WBL will 
break down if the Fermi energy is very close to a band 
edge. 

Within this model of a gold chain, it is probably most 
physical to occupy the atoms with a single electron per 
site. However, in a more realistic description the s and 
Pz orbitals hybridize also with lower-energy d orbitals. 
Also, the real systems are not translationally invariant, 
and charge neutrality need not be fully satisfied locally. 
Such things can complicate the picture enormously. Fur- 
thermore, the arbitrariness of the charge-neutrality pro- 
cedure in the full spd model below renders the position 
of the Fermi energy with respect to the electronic struc- 
ture of the wire somewhat uncertain. Thus, as we shall 
see, all of the cases shown in Fig. |S1 are actually good 
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FIG. 5: The top panels show the conductance vs. voltage 
G{V) for the spz chain model with parameters as in Table 
and the lower panels show the derivative dG/dV. Here the 
curves have been shifted by integer multiples of 0.01 or 0.005 
for clarity. The left panels are at T = 1.0 K, while the right 
panels are at T = 4.2 K. The vibrational modes are as in Fig. 
1 of Ref.lM 



descriptions of our results for the full spd model. 

Although at low temperature there can be several 
peaks, at higher temperatures the steps fuse together, 
since the vibrational frequencies of the different modes 
are quite close to each other. The fusing effect of the 
steps would be further enhanced, if we introduced a larger 
broadening rj for the vibrational modes - here t] = 0.002 
meV <§; ksT. The small rj also affects the results in 
another way: since the "lead coupling" of vibrations is 
small, the wire heats, and the distribution functions Na 
differ strongly from the Bose distribution. The signature 
of this heating is the steep downward slope of the G{V) 
curves at high voltage [l9l |. 

We have also tested the effect of n'th-nearest-neighbor 
hoppings and addition of nonorthogonality, but these 
have no essential qualitative effect on the results. These 
will be taken into account in the full spd parametrization, 
which we now turn to. 



VI. FULL spd PARAMETRIZATION 

In this section we describe a more realistic tight- 
binding approach to the problem. We use the nine- 
orbital spd parametrization of Papaconstantopoulos et 
al. miil il El. This type of spd TB approach is 
known to reproduce very well some nontrivial ab initio 
results, like the numbers of conduction channels and the 
formation of zigzag Au chains ^38, j|9j . Thus we can be 
confident that the method gives at least good order-of- 
magnitude estimates for all of the quantities which we 
shall be interested in. 

However, since the parameters are extracted from first- 
principles bulk calculations, they cannot be exactly cor- 
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FIG. 6: Dimensions of the unoptimized geometry with N^h ~ 
4. Here a = 4.08 A, the fee lattiee eonstant. Only the eoor- 
dinates of the Nch ehain atoms are optimized. 



rect for atomic point contacts, where the important 
atoms of the structure are significantly less coordinated 
than in bulk. It has thus become customary in the 
method to "correct" the parameters in the central cluster 
in order to satisfy local charge neutrality. Doing this typ- 
ically brings the central cluster levels better in resonance 
with the lead orbitals. We compute the charge with the 
so-called MuUiken population analysis (see Appendix IB|| . 
and only shift the on-site energies of the Hamiltonian. 
Tests with other ways of achieving neutrality give very 
similar results 28]. Furthermore, we have compared to 
results obtained without charge neutrality, and again find 
that there is not much qualitative difference, although 
charge neutrality yields conductances closer to Go on 
average. Also, without charge neutrality there appears 
to be a tendency for seeing lower-frequency vibrational 
modes in the conductance curves. 

The results of this paper are generated using finite lead 
broadenings Jl,r ~ 1-0 eV. The limit 7L,i? O"*" can in 
principle be taken without affecting the results in any 
essential way. 



A. Geometry optimization and vibrational modes 

We consider two types of ideal geometries, the "A" 
and "B" ones, shown in Figs. ^ and |21 respectively. As 
mentioned, the leads are assumed to be of fee type with 
the [001] axis in the transport (or z) direction. Before 
geometry optimization, the chain atoms are positioned 
as described in Fig. El The "length of the wire" Lch is 
defined as Lch = a-\-d{Nch — 1), where Nch is the number 
of atoms in the chain, d the distance between them, and 
a = 4.08 A is the equilibrium lattice constant of the 
bulk fee lattice. We only optimize the geometry of the 
Nch chain atoms — also in geometry B which has the 
"pyramids". Thus, although the interatomic distances 
change slightly from those of Fig. 1^1 Lch remains fixed. 

To estimate the total energy E(Rk), we simply take a 
cluster which includes the wire and some atoms from the 
leads, solve for the electronic eigenstates £„, and then 
occupy the states according to charge neutrality. This 
energy, as a function of the 3Nch wire coordinates, is 
then optimized with standard library routines. As known 
previously ^8], it is often energetically favorable for the 
gold chains to exist in a zigzag-like pattern instead of a 
linear one. Only after a sufficient amount of stretching 



(i.e., with a larger d) does the linear configuration become 
stable, after which it remains linear until the wire breaks. 
We find that the maximum d at breaking is, depending 
on the geometry, typically something between 2.70-2.85 
A. There is no clear trend with increasing Nch- 

After the geometry is optimized, the energy function is 
used to compute the Hessian matrix Ti.. The eigenvalues 
ka {a — 1, ... , 3Nch) are all positive, and the vibrational 
frequencies are simply given by uja — \/ka /M. With 
both geometries, A and B, we obtain quite similar vi- 
brational frequencies and modes. For a linear wire, the 
modes can be classified as longitudinal or transverse in 
character. The highest-frequency modes are then always 
longitudinal ones, and the highest of them is of the ABL 
type M- 



B. Elastic transmission 

Perhaps the most characteristic experimental property 
of gold chains is that they appear to have a conductance 
very close to the quantum of conductance Go. Thus we 
shall briefly comment on the elastic transmission proper- 
ties of the chains in our calculations. The present TB 
method was previously only used with bulk distances 
d w 2.885 A between all the atoms |3- In this case, 
we find very similar transmission functions To(e) for our 
geometry B. These are often characterized by very long 
plateaus (of widths up to 2 — 3 eV) around the Fermi 
energy ep, where 0.95 < To ^ 1-0. The same is true for 
the results with and without charge neutrality. 

However, when the geometry is optimized, the wide 
transmission plateaus close to one are replaced by larger 
oscillations. Still, at the Fermi energy, there is usually 
only a single open channel, which consists of s, Pz, and 
(i322„r2 orbitals. Sometimes, a small contribution is seen 
arising from a second channel, involvin g th e other p and 
d orbitals, as will be discussed below [23]. The trans- 
mission around ep varies between 0.7 ^ To < 1.0. The 
present method is known to reproduce experimental con- 
ductance histograms rather well 40]. In particular, the 
conductance peak somewhat below G/Gq = 1 is a very 
robust feature. 

We determine by charge neutrality in the leads (or 
bulk). Tests with shifting its value from this position by 
^ 0.5 eV (which simulates a gate effect) showed no sig- 
nificant qualitative effects on the results. The position 
of the Fermi energy with respect to the local electronic 
structure is still very important. This is because, be- 
sides the elastic transmission, the Fermi energy also fixes 
the "Fermi wavevector" , which again determines what 
vibrational modes can be excited. In the present TB 
method, the use of bulk parameters and the charge neu- 
trality procedure introduces some uncertainty in relation 
to this point. 
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FIG. 7: A sketch of the vibrational eigenmodes of a linear 
four-atom gold wire, with energies corresponding to d = 2.62 
A in geometry A. The longitudinal modes (L) are all non- 
degenerate, whereas the transverse modes (T) are all doubly 
degenerate. 



C. Longitudinal and transverse modes 

Let us first discuss the basic observations using a sim- 
ple example, namely, a linear chain of four atoms in 
geometry A. A schematic illustration of the vibrational 
modes is shown in Fig. □ for d = 2.62 A. There are four 
longitudinal modes and eight transverse modes. How- 
ever, due to the fourfold rotation symmetry of the geom- 
etry around the axis of the wire, the transverse modes are 
all doubly degenerate. The zero-bias conductance is due 
to two partially open channels. The main contribution 
(about 98% of Go) is due to a channel (Ci ) formed from s, 
Pz, and (i3^2_^2 orbitals, which have the symmetry of the 
geometry. In addition, there is a small (less than 1% of 
Go) contribution from a second, doubly degenerate chan- 
nel (G2), which consists of dxz, dyz, Px, and Py orbitals, 
which have a lower symmetry. Thanks to the symmetry 
of the Gi channel, only longitudinal modes have a finite 
coupling constant in its subspace (A^^ (^J. In the sub- 
space of the G2 channel, also the transverse modes have a 
finite coupling {^02,02) ■ Thus we might expect that also 
the transverse modes give a small signal in the current. 

Figurc|Slshows an analysis of the contribution from the 
different modes to the differential conductance G{V) — 
dl/dV. We divide this conductance into three parts 
according to the three current contributions: G{V) = 
GoTo + SG^nei{V) + SGec{V). Here Go = 2eyh, 6G,nei = 
dlinei/dV, and SGec = dSIec/dV. It is seen that SGinei 
gives always positive contributions to the conductance 
steps, while SGec gives negative ones. As expected, we 
find that there is a finite step in both SGmei and SG^c 
due to all of the vibrational modes, also the transverse 
ones, although the latter are quite small. However, the 
contributions of SGi„ei and SGec for the transverse modes 
almost perfectly cancel each other, such that only steps 
due to the longitudinal modes are seen in the total G{V). 



FIG. 8: Decomposition of the conductance G{V) into GoTo, 
an "inelastic" contribution ddnei , and an "elastic correction" 
SGec for a four-atom wire. The geometry and the labels (a)- 
(h) correspond to those of Fig. [7| Other parameters are T — 
0.01 K, r; = 0.002 meV. The solid step-like curve shows SGinei, 
and the dashed one shows SGec- The increases of SGi„ei due to 
transverse modes are exactly canceled by decreases in SGec- 
The inset shows the elastic transmission (dashed line), and 
the total conductance G{V) (solid line). In G{V) only drops 
due to longitudinal modes are seen. 

This cancellation is apparently due to the exact fourfold 
rotation symmetry, and the mirror symmetry with re- 
spect to the plane cutting the wire in the middle. In less 
symmetric geometries the transverse modes can also give 
finite contributions to G{V)- 

In the case of zigzag wires, the distinction between lon- 
gitudinal and transverse modes does not really exist, and 
all modes are always seen as steps in G{V)- An exam- 
ple of this is shown below. Furthermore, if the elastic 
transmission To is very small, then also the transport co- 
efficients r^'^ and T^'^^^ tend to be small, since they all 

~ r ~ a 

depend on the matrix G TjiG T^. In this way, for exam- 
ple, it may also be possible to have large positive steps 
in G{V), but we never see them for the charge-neutral 
gold wires. For other materials, the situation may be 
different. 

Thus, we find that the conductance features depend in 
an intricate way on the symmetries of the geometry, the 
symmetries of the vibrational modes, the coupling con- 
stants, as well as the symmetries of the electronic states 
which are relevant at the Fermi energy. 

D. Conductance curves of linear wires 

Here we discuss in more detail, how our conductance- 
voltage curves for linear wires look like. Figure |51 shows 
an example for geometry A with a wire of Nch — 4 atoms 
while Fig.^]is for a wire of Nch = 11 atoms in geometry 
B. In both figures, the left-hand panels are calculated 
at r = 4.2 K and with a small rj, such that they are 
more or less in the externally undamped regime. The 
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FIG. 9: Comparison between theory and experiment for 
Nch = 4 in geometry A. The sohd and dashed curves cor- 
respond to theoretical results for d — 2.54 A and d = 2.68 A, 
respectively. On the left-hand panels T = 4.2 K, whereas on 
the right-hand panels these curves have been broadened with 
a larger temperature T = 12 K; in both cases rj = 0.02 meV. 
The experimental results LI and L4 ( x jcorrespond to 
the notation and results of Fig. 1(d) of Ref. [13 with V > 0. 
They are obtained for a 7-atom chain at T = 4.2 K. 
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FIG. 10: Comparison between theory and experiment for 
Nch = 11 in geometry B. AU resuhs are at T = 4.2 K. The 
solid and dashed curves correspond to theoretical results for 
d = 2.64 A and d = 2.78 A, respectively. On the left-hand 
panels 77 = 0.002 meV, whereas on the right-hand panels the 
curves have been broadened with a bath-coupling r; = 5.0 
meV. The experimental results LI and L4 (x) are as in 
Fig. El 



right-hand panels show two examples of the experimen- 
tal results for a wire of approximately 7 atoms taken at 
the temperatm'e T — 4.2 K J_^J. Comparing these to the 
theoretical curves on the left-hand side, one immediately 
notices that if the conductance drop is to be due to a 
single mode, then the ~ 5 meV width of the peak in the 
experimental dG/dV cannot be explained by tempera- 
ture alone 19]. On the other hand, the energy distance 
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FIG. 11; Zero-bias conductances, voltage positions V^h, and 
the local maxima of IdG/dVI for different wire lengths L^h 
and atom numbers Nch, indicated by the numbers. The 
results are obtained for linear wires in geometry B, with 
rj — 0.02 meV. The dashed lines show the limits beyond which 
the wire breaks upon further stretching. On the left side of 
each curve, the wire has a non-linear form, and the peak struc- 
ture may be different (see text). The solid lines correspond 
to the main peaks and the dotted lines to the second-largest 
ones. Thickening of the lines indicates that the peak consists 
of two close-lying modes. 



between the vibrational modes is rather large > fc^T, 
so that at T = 4.2 K, a peak consisting of several sub- 
peaks can usually be easily recognized. For example, the 
highest-frequency peak in Fig. IIUI actually consists of two 
peaks, and it is still not wide enough. Also in Fig.|51there 
are at least three separate peaks visible. 

Thus we conclude that in the experiment there are 
probably other broadening mechanisms at play besides 
temperature. In the right panels of Fig. |51 we compare 
the experiment with a theoretical result broadened by a 
higher temperature, while in the right panels of Fig. IIUI 
we use the parameter rj to broaden the peaks. In the 
latter, the system is already in the externally damped 
regime, with very little local heating: in addition to the 
broadening, the damping is signified by a smaller slope 
after the drop. In either case, the peaks due to individual 
modes are smoothed out to form a single one, with a 
width comparable to that seen in experiments. In this 
way, it is possible to obtain a rather good quantitative 
correspondence between theory and experiment. 

We have also studied systematically how the main fea- 
tures of the conductance- voltage curves vary when linear 
chains with atom numbers 3 < Nch < 11 are stretched. 
The results for T — 4.2 K and a small 77 are plotted in Fig. 
im Here we show the zero-bias conductance G{V = 0), 
the voltage positions Vph of the main peaks in dG/dV, 
and the corresponding local maximum values of |dG/dT^| . 
The peak height is not the same as the peak area which is 
plotted in Ref.[ll but is roughly proportional to it, since 
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the width of the peaks is fcsT, where the temperature 
is T = 4.2 K. 

In this figure, many features can be recognized. The 
G{V = 0) values fluctuate between 0.85 and 1.0, but 
there is no clear "parity effect" — or at least the effect ap- 
pears to reverse its direction after Nch = 8. The present 
method most likely does not describe such parity effects 
correctly. However, also different ab initio approaches 
are known to give conflicting results 0- In some calcu- 
lations, the transmission has been found to oscillate also 
with the stretching of a wire with fixed N^h US • 

The positions of the peaks move to lower voltages when 
the wire is stretched, as expected from the "softening" of 
the bonds, and the resulting decrease in the vibrational 
frequencies. There is also a clear trend toward lower fre- 
quencies with increasing Nch- These findings are similar 
to what is seen in the experiment |l^ . Also similarly, 
the peak heights increase with stretching, and with in- 
creasing Nch, although the increase with wire length Lch 
is not obviously linear as for the simple chain models 
of Sec. El 20] . There is also a correlation between the 
zero-bias conductances and the peak heights: when the 
conductance is low, also the stretching behavior is rather 
anomalous (in particular in the cases of 5 and 10 atoms.) 

The most visible difference between the results of Fig. 
111! and the experiment is that we consistently see sig- 
natures of several vibrational modes: typically there are 
two peaks visible in dG/dV. However, the higher peak is 
always at a larger voltage and, as the number of atoms 
Nch grows, the secondary peaks become less and less sig- 
nificant. For example, for 11 atoms there is essentially 
only a single peak visible (see Fig. I10|) . This, however, 
is due to two close-lying modes: the highest-frequency 
"ABL" mode, and the one next to it in frequency. As 
explained above, this discrepancy of several peaks can be 
corrected by increasing the parameter rj. Note also that 
this behavior was already present in the chain model of 
Sec. El 

We also see that the largest conductance drops are sys- 
tematically at too high voltages compared to experimen- 
tal values Vph ~ 10 — 20 mV for a 7-atom chain. This 
is not surprising, given the simplified way in which we 
compute the vibrational modes. The frequencies of the 
vibrational modes might be lowered, if we also allowed for 
the motion of atoms outside of the chain. In other words, 
it is possible that the lead coupling, done in a proper way, 
would lead to a "redshift" of the frequencies. As noted 
above, the electron-vibration coupling gives such a red- 
shift [23 , but this effect may be too small to explain the 
discrepancy. 

Although the steps in conductance are almost always 
downward when G{V — 0) is close to Go, sometimes 
also weak increases in the conductance at low voltage can 
be seen. These appear to be related to the longitudinal 
"center-of-mass" mode. For linear chains, we do not find 
any significant contribution from transverse modes, as 
explained previously. 
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FIG. 12: The difference between zigzag and linear wires of 6 
atoms in geometry A at T = 4.2 K and with rj = 5.0 meV. 
The numbers in the legend indicate the parameter d/A, which 
describes the amount of stretching. In the left panels the wire 
has a zigzag character, while in the right panels it is linear. In 
case of the zigzag-wire, practically all of the vibrational modes 
contribute to the observed steps, but they form a clear double- 
step structure. In the linear wire, the elastic transmission 
varies a lot under stretching in this example. 



E. Zigzag wires 

We have also studied briefly the zigzag chains shown 
in Fig. n Figure 1121 shows a comparison between the 
signatures of a zigzag chain and a linear chain in the 
conductance- voltage curves G{V). In this example, a 
chain of Nch = 6 atoms in geometry A, and a large broad- 
ening of the vibrational modes {rj = 5.0 meV) was used. 
However, very similar results were found for various dif- 
ferent atom numbers. For small values of d, the wire has 
a zigzag character, but at d « 2.42 A, the wire becomes 
linear. 

For the zigzag chain, there are two well-separated (se- 
ries of) conductance drops. The one at low bias is higher 
for small d's. As the chain is stretched, the height of the 
lower-frequency step decreases, while that of the higher- 
frequency one increases. For the linear chain, there is 
essentially only one high-voltage drop, which is actually 
due to two different longitudinal modes. In the case of 
a zigzag wire, the vibrational modes have complicated 
symmetries, and practically all of them are contribut- 
ing to the current steps. However, here their identities 
are completely smeared out due to the broadening. This 
calculation provides a clear prediction of how the zigzag- 
to-linear transition may be seen in experiments. 



VII. CONCLUSIONS AND DISCUSSION 

We have studied the onset of dissipation by excitation 
of vibrational modes in atomic gold wires, using a tight- 
binding model. We have studied in a systematic way 
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how the stretching of wires with different atomic num- 
bers affects the conductance steps, and find a reasonable 
agreement with experiments. Previously such a study 
has only been done within simple chain models p3(- We 
have also considered two different geometries, which yield 
qualitatively similar results. Finally, we studied the con- 
ductance signatures of zigzag-type wires, and predict a 
double-step structure in contrast to the single step of lin- 
ear chains. 

Our results for the linear chains agree rather well with 
experiments and previous ab initio calculations, apart 
from the incomplete "mode selectivity" . In this context, 
we have pointed out the importance of taking into ac- 
count the broadening of vibrational modes due to their 
coupling to the leads, especially in the limit where the 
vibrational mode distribution is assumed to be strongly 
damped. We derived equations in the wide-band limit, 
which take this into account in a phenomenological man- 
ner. The wide-band limit combined with the lowest-order 
perturbation approach (presented in Sec. IIVI and Ap- 
pendix appears to provide a sufficiently good descrip- 
tion of the phenomenology of electron-vibration interac- 
tion in atomic wires. To make further progress, more de- 
tailed calculations of the lead-coupling of the vibrational 
modes are needed. 

As explained above, the condition which determines 
what modes yield conductance drops is essentially that 
of momentum conservation. We find numerically that 
the momentum-conservation idea works well also for fi- 
nite wires with 10 atoms or more. However, in shorter 
wires, the electronic structure is very complicated, and 
the connection between symmetries of excited modes and 
those of the electronic states is hard to analyze. 

Even if there is approximate momentum conservation, 
there is no fundamental reason why only a single mode 
would be excited, or that the highest-frequency (or other 
ABL) modes should necessarily be involved. This re- 
mains true also for very long wires, since at the same 
time when the momentum conservation becomes more 
and more accurately satisfied, also the energy density 
(and hence momentum-space density) of the modes in- 
creases. It is thus more likely that a large "wave packet" 
of several nearby phonon modes is always excited. How- 
ever, the charge neutrality of the wire indeed appears 
to be favoring the highest-frequency modes, just as pre- 
dicted by the simple chain models discussed in Sec. 
Nevertheless, it is important to note that small differ- 
ences in details of the geometry or the implementation 
may affect numerical values of results significantly. The 
band structure of an infinite, linear gold wire is already 
quite complicated [2^ [sH , such that a small error in the 
Fermi energy, and hence of the Fermi wavevector, will 
immediately shift the resonance to slightly different vi- 
brational modes. 

Studying different materials (Ft and Ir) with this same 
method would be straightforward in principle. However, 
it seems that the parameters available for these materi- 
als are not very good for geometry optimizations of the 



wires. This is because the overlap matrices easily lose 
their positive definiteness when the validity range of the 
parametrization is exceeded (with the gold parameters 
|25], such problems never appeared). Thus, as already 
implied in Ref. 28, one should probably use more gen- 
eral ab initio methods, or at least parameters which have 
been specifically fit to work for chain geometries with a 
large span of interatomic distances. Compared to ab ini- 
tio methods, TB still has the clear advantage of being 
computationally efficient. 
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APPENDIX A: CALCULATION OF THE 
COUPLING MATRIX ELEMENTS 

Here we describe our method of computing the ma- 
trix elements M^J^ = (i|Vfc^iJ|j), where the derivative is 
with respect to the components /i = x,y,z of the ionic 
coordinates Rk- In an implementation making use of 
TB parametrization, one has no direct access to either 
the basis states or \i), or the Hamiltonian H — only 
the representation Hij = {i\H\j) and the overlap matrix 
Sij = {i\j) as a function of Rk are known. Thus, the ma- 
trix elements must be calculated more indirectly. Below 
we sketch a derivation, which follows, in some sense, the 
ideas of Ref. ^3- The derivation is not exact, as we only 
consider an isolated central cluster. 

Let \i),\j),... denote the atomic orbital (AO) basis 
states, and . . . the electronic eigenstates (molec- 

ular orbitals, MO) of the central-cluster electronic Hamil- 
tonian H. The eigenstates satisfy 



{a\f3) = Sap. 

If the expansion of the MO's in the AO basis is denoted 
a) — J2i \i)Cia, where Cia — {ca)i then we have the 
matrix equations 

Hcq, Sc^Cct 

Let us also note the form of the completeness relation of 
the MO's J2a = 1 in the AO basis: EaCacj, = 

while for the AO's themselves X^^- l«)(5'"^)<j 01 = 
1. Since the basis states also move with the ions, we 
should write more carefully \i{Q)), \a{Q)), H{Q), Hij{Q) 
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etc., where Q is a shorthand for the ionic displacements 
Qk = Rk — -Rfe'''- The point in the electron- vibration 
coupling is that moving an ion will induce a perturbation 
of the form 



H 

S 



H 



-H'Q 
S'Q 



(A3) 



This is the final expression for the matrix elements. The 
overlap corrections due to the nonorthogonal basis often 
turn out to be rather small in practice, an the orthogo- 
nal result M — H' is a reasonable first approximation. 
The presence of the corrections tends to make the con- 
ductance steps slightly larger. 



and one can calculate its effect on the state vectors 
and eigenenergies Ca by means of simple first-order per- 
turbation theory. Here " ' " denotes a derivative with 
respect to Q. 

Let us expand the basis states as \i{Q)) ~ \i) + \i')Q, 
where \i) now denote the unperturbed basis. The matrix 
element which we are looking for is really the quantity 



My- = {i\H'{Q)\j)\ 



(A4) 



This may be obtained by considering the expansion of 
{i\H{Q)\j) in Q, since 

{imQ)\j) - ^..(0) + m'mj)\Q^^Q. (as) 

By inserting 1 ~ \a[Q)){a{Q)\ one easily finds 

mmj) = ^(*|a(g))ea(g)(a(Q)|j). (A6) 



Then, by inserting \a{Q)) = Y.i \i{Q))Cia{Q), expanding 
ea{Q) « + e'^Q, CiaiQ) « Cia + C-^Q, and equating 
terms linear in Q 



kl 



' '-^ik '-^ka^a'-^la'-'lj ~r >~'ik'~^kata'~^la^lj 



(A7) 



Here we introduced the one-sided overlap derivatives 
[S'^%j = [S'^%j = (%•'), which satisfy S' = 

+ and = [S'(i)]t 41J. 

The result of Eq. ^Mil may be simplified considerably 
by inserting the expressions for and C^'^,, which are 
easily derived. First, if we denote 



V(l,2) 



kl 

'-ka>=>ij 

kl 



kl 



and so on, one may show that 



c'(l)^ 



'(2) 
Q/3 • 



Then, using the completeness 

Q fc C'fcQ C*^ = J2aJ2k^ia ^ka ^kj 

may transform back to the AO basis: 



EE 

a kl 



{ 



c'(l)/^ ^ Q ] Q ^ C^2) \ 



V(2) 



10, 



(A8) 

(A9) 

relations 
6ij one 

(AlO) 



APPENDIX B: POPULATION ANALYSIS AND 
DERIVATION OF THE CURRENT FORMULA 

Although the expression for the current is well known, 
its derivation in the nonorthogonal basis is not entirely 
trivial 42j . The current is the time derivative of charge 
transported from L to C and on to R. However, while 
the total charge of the full system is well defined, the 
partial charges of the distinct regions are not — instead, 
there arc different ways of doing "population analysis" . 
To be more self-contained, we present here a compact 
discussion of these issues. 

Let us consider, for simplicity, a single orbital per 
atomic site i. The total charge (particle number, ac- 
tually) is given by 



kj. 



3,k 



(Bl) 



where Q ~ ^ d\ Skj dj and we define the density ma- 
trix Pjk = {d\.dj). Partial charges may be introduced for 
exarnple with the following "MuUiken" population anal- 
ysis [43 



Q = Ql + Qc + Qr^ J2 P^kSkj 
jeL,k 

+ E PjkSkj + X! PjkSkj- 
jeC,k j£R-k 

One may also define the quantities 



(B2) 



- (Qn) - E PjkSkj, 



j.ken 



(B3) 



where = J2j,ken'H:^k]dj, and CI L,C\R,C + 
R,C + L. Here C + L {C + R) refers to the combined sys- 
tem involving C and L (C and R) regions. The charges 
Q'l c r ^'^^ good approximations to Ql,c,r if the regions 
are all large, since the corrections are proportional to the 
surface area of the interfaces. 

What we want is the particle current J through a 
boundary between two regions of space, L and C, for 
example. We expect it to satisfy 



2J 

Let us rewrite Q'^ - 



d_ 

dt 



{Q'l 



IC+R 



(B4) 



Nn 



Y.j.k&n,j^kdlSkjdj, where 



and N, 



d\dj for 



L,C + R. 
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Now, in the Heisenberg picture the dj operators sat- 
isfy the equations of motion ihdj = [dj , H^] and hence 
ih£i 'l2k'^jkdk = Hjkdu- Using ideas similar to Ref. 



144 we find 



dt 



Next, considering the following quantity, we notice that 
all hopping contributions except Hlc and Hlr cancel: 



E 



jeL,kec+B. 



dj— (.Hjk - Sjkih-^ ) dk + h.c. 



E 



jeC+R,keL 



dl^[H,k~S,kih^^jdk+ h.c. ^gg^ 



d 



E 4 Sjkg-^dk+ h.c 



j,keL,k^j 



E dlSjk-g-^dk+ h.c. 

j,keC+R.,k^j 



When the last two overlap terms are moved to the left- 
hand side, it becomes exactly what we called 2J. Thus 
the expectation value is 

2J = 2(J) 

= E UH,k~S,kih^){d]it')dk{t)) 



jeL,kec+B. 



E 



jeC+R,k£L 

c.c. 



dt 

,^ H,k-S,k^h^^ 



{d]{t')dk{t)) 



(B7) 



Now, assuming that H j^j^ = Slr = and using 
i{d]{t')dk{t)) = G+-{t,t'), we finally have 



APPENDIX C: NEGF FORMALISM: 
TECHNICAL DETAILS 

Our notation differs slightly from what is the standard. 
In particular, our functions are equal to the 

functions of Refs. 19, 32], for example. Let us write the 
definitions for the most important electron propagators 

G:^{t,t') = -i{{d.,{t),d]{t')})0it -t') 

G1^{t,t') = i{{d,{t),d]{t')})0{t' -t) 



Gt-{t,t')^i{d]{t')dm 

G^+{t,t') = -i{d,{t)d]{t')). 



(CI) 



Similar expressions hold for phonons, but with the re- 
placement {•, •} — > [•,•]. Below, all Green functions ap- 
pear Fourier transformed with respect to t — t' into an 
energy representation. 



1. Electron propagators and self-energies 

All expressions for the relevant electron Green func- 
tions follow from the "Dyson" and "Kcldysh" equations 

G^^ = (1 + G'-E'')g^^{l + S"G") - G^'S^^G"', 

(C2) 

where the upper or lower signs can be chosen. Here the 
g denote the Green functions for an uncoupled central 
part in the absence of electron- vibration interactions and 
S*"'*^ are the sums of all self-energies containing the 
effects of both (see below) . The uncoupled functions are 



9+- = -fci9'' -an 

9-+ = -{.fc-m9'-9n- 



(C3) 



J(t) = --RcTr 

h 



HLc~SLcih-]G+l{t,t') 



(^HcL-ScLih^^ Gtcit,t') 



(B8) 



The charge current is then obtained as / = — 2e J, where 
the factor 2 is the spin degeneracy. An analogous deriva- 
tion may be carried out for the current over the C-R 
boundary. 

In stationary state all the propagators only depend on 
the time difference t — t' , and this result may be Fourier 
transformed into an energy representation, where one has 
to replace ihd/dt — > e. 



Here is the equilibrium Fermi distribution. Note that 
7c = i[(9'°)^^ — (9'^)^^]) which is an infinitesimal quan- 
tity. 

The functions where the lead coupling is taken into 
account but electron-vibration coupling is still neglected 
are given by 

G = [eScc - Hcc - ^ ^fi] ^ 
G^ = -G'-(S+- + S+- - ncfc)G'' (C4) 
G^ = + - i7c(/c - 

and G — {G Y . Here the lead self-energies and lead 
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Green functions for the L side are given by 

glL^[{e + hL/2)SLL-HLL]-' 
Tl = i(Sl - S2) (C5) 

= tcLgth^Lc = -iri/i 

= tCLQlttLC = -irL(/L - 1), 

where tie = H lc — ^S lc and so on, with similar expres- 
sions for R side. The infinitesimal ')q is only needed for 
recovering the correct result in the limit where the self- 
energies are taken to zero — it may be neglected here. 
The parameters ^l,r are positive infinitesimals, which, 
however, can be used to introduce a finite broadening of 
the lead eigenstates. 

These are enough for calculating the elastic current 
in the absence of electron-vibration coupling. The full 
Green functions including the effects of this coupling are 



G 



-+ 



- ndfc - ^)]G'' 

(C6) 



and G" = (G^)^ To second order in the coupling con- 
stant A", the electron-phonon self-energies are 



e—vib 



''e—vib \ 



a 

+ i?j:K)[A"GT±(e-L.i)A"] 
-A"Tr[G+-K)A"]i?;(0)} 

-f Z?SK)[A"GT±(e-c.i)A"] 
-A"Tr[G+-K)A"]i?^,(0)}, 



(C7) 



where the upper or lower signs are chosen. 



2. Phonon propagators and self-energies 

By our definition of the phonon propagators, the un- 
perturbed ones (those in the absence of a lead coupling 
and electron- vibration coupling) are given by 



1 



1 



e-eQ-|-i7//2 e + ea+i'q/2 
2e„ 



— -|- irje — rf /A 
d+"(e) = -~2iT'm{e)pa{e) 
d-+(e) = -2^i[n(e) + l]p„(e). 



(C8) 



Here = fvjJa are the bare vibrational energies, n(e) 
is the Bose distribution function. The quantity pa = 
— Imd^/TT is the bare phonon density of states (DOS), 
given by Eq. ©. This becomes Pa{f) = '5(e — £») —5{e + 
Eq), as rj 0+. Note that the Green functions satisfy 
rj = i[(d'')-i - (d'-)-i]. 

Now, to be symmetric with the discussion for the elec- 
tron propagators, the next step should be the introduc- 
tion propagators "D" which contain lead self-energies. 
For the proper calculation of the lead self-energies, we 
should probably change to an atomic-displacement ba- 
sis. Instead of doing this, we model the lead coupling by 
giving finite values to the infinitesimal quantity rj, which 
will broaden the phonon density of states ■ 

The full phonon propagators appearing in the electron- 
phonon self-energies are obtained from the "Dyson" and 
"Keldysh" equations for phonons 



(C9) 



For self-energies, or "polarizations" , we use the second- 
order approximations 

n^^(e)=i J ^Tr[A"G±^(.;OA"G^±(cc-i-e)] 

n^W = -i/ ^Tr[A"G±^(.;i)A"G''(c.i -6) 
A"G'^(a;i)A"G±=F(t^^ 

noW = -i / ^Tr[A"G±^(c.i)A"G^(c^i-6) 
+ A"G'^(c^i)A"G±T(c^i -e)], 

(CIO) 

where we have dropped some unimportant zero- frequency 
terms. Again either the upper or the lower signs must be 
chosen. Note that are purely imaginary and satisfy 
the symmetry 11+^ (—e) = n^+(e). 

Equations IjClOp close the system of equations, and 
we are done. However, from a physical point of view, it 
is interesting to develop the equations slightly further. 
Using the symmetry — _D° = D~~^ — -D^~, one finds 
that Eqs. ljC9p may be rewritten in the form 



2e„ 



£2-4 +i,ye- 7^2/4 -2e„H-(e) 



D+-{e)^~2TriN^{e)paie) 
D-+{e) = -2n\{N^(e) + l)p^{e), 

where we define the phonon density of states 
p„(e) = -ilmi?^(e). 



(Cll) 



(C12) 



which satisfies /Oa(— e) = — ^^(e). The quantity Na{e) 
is the energy distribution function of the vibrational 
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quanta. In equilibrium Na{e) = n{e), the Bose distri- 
bution. In general Na{e) = n{€) + 6Na{€), where SNa{e) 
is a voltage-dependent non-equilibrium correction. 
By reshuffling the Keldysh equations one may write 



D-+{e) = -i?^[i(l + n)er^/e^ + Il-+]D^^. 



(C13) 



Since jf^P = Im£'J^/(Ininj^(e) — rye/2ea), comparing 
these with Eqs. (|("1H) it is easy to obtain explicit expres- 
sions for the distribution function N^: 



llmn+-(e) H-n(e)rye/e^ 
2 Imn-(e)-r;e/2e„ 
llmn-+(e) + [l + n(e)]?7e/e„ 
2' 



(C14) 



1. 



Imn^(e) -77e/2ea 
To get to the last line, we used IT^ — 11° ~ ^(B^ 



nj~). Note that we have not made any approximations 
to get to this result. Using Eqs. HC10() . one finds that the 
following symmetries are valid: Na{—e) = —[Na{e-) + 1] 
and SNa{-e) = ~SNa{e). 



APPENDIX D: CURRENT CONSERVATION 

The inelastic parts of the current in Eqs. Q look very 
asymmetric in their L or R indices. Nevertheless, in sta- 
tionary state the currents calculated at L and R bound- 
aries should be equal: — . There also appears to 
be some confusion as to what sort of approximations are 
needed for current conservation |,19„ .36] . Here we outline 
a proof of this property for our approximation. For I^i it 
is obvious, so we only consider the inelastic current. 

Inserting the self-energies from Eq. l|C7p into the ex- 
pression for /j^jgj [Eq. 0], one finds 



inel 



E 



2e 
Tr[G'^(e 



de f duji 
2^ J ~2 



— 2Trpa{LJi)\ 
TT L 



^1, 



^)A"G''(6 + y)ri(6- 



UJl 



X [fde + ^){N^{uJi) + 1) - {fde + ^) + iV„(c^i))/L(^ " y)] 



-Tr[G'^(e-^)A"G''(6 + ^)ra6- 



X -I- y )(A^„(a;i) + 1) - (/l(6 +'^) + N^{cji))fn{e - ^)]}. 



(Dl) 



Here the first term may be shown to vanish as fol- 
lows. Since Fermi and Bose functions / and n satisfy 
f{x)[n{y) + 1]- [fix) + n{y)]f{x - y) = 0, we have 



fL{x)[N^{y) + 1] - ihix) + No,{y)]fL{x - y) 
^SNMhix)- fdx-y)]. 



Now, noting that 5Na{—e) — —6Na{e), and e) = 
~Pa{e), and additionally assuming that all matrices un- 
der the trace are symmetric (G'' = C, [F^]^ = Tl 
etc.) it is seen that the integrand is odd and thus the 
energy integral vanishes. This is true also if (57V„(e) has 
a 1/e divergence at e = 0, since the product in Eq. 
(|D2|I remains finite. 

The second term in Eq. (jPip is clearly symmetric upon 
interchanging L and _R, and changing the overall sign. 
Note that the proof does not rely on things like a mirror 
symmetry of the geometry, and remains unchanged if we 
replace G^ —>■ G . Thus our expressions are always "cur- 
rent conserving", as defined by the condition 1^ = 1^. 



APPENDIX E: GENERAL PERTURBATIVE 
CURRENT FORMULAS 

Here we shall write down the perturbative current for- 
mulas of Eqs. (O in a slightly different form, which will 
make the so-called wide-band approximation more trans- 
parent. However, here we shall be making no approxima- 
tions in addition to the second-order perturbation theory 
which we have already introduced. Nevertheless, it must 
be stressed that, due to the second-order approximation, 
all higher-order terms in the following expressions are 
strictly speaking not warranted. This concerns the ap- 
proximations made for the vibrational density of states 
Pa and the distribution Na, which should in principle 
both be of zeroth order in the electron-vibration cou- 
pling constant. In the approximations which we use for 
Na, this is not necessarily the case. Thus our approach 
is not purely lowest-oder perturbation theory. However, 
the self-consistent Born scheme (SCBA) is not really any 
better is this respect. 

The inelastic current /i„, and the elastic parts I^^ 
and 6Iei are obtained from Eqs. Inserting the self- 
energies [Eqs. (|C7|) ] into these formulas, they may be 
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rewritten as follows: 



/° =^ I deToie)[hie)-fRie)] (El) 
6Iei = - / d6^ " / dc.ip„(.^i)|T-(e,c.i)7V„(ac.i) +T-^(6,c.i)/i(6,„) + r-«(e,c.i)/fl(e.„)| 

Li •'0 



a o"— ±1 



(E2) 

X [/L(e) - fR{e)] + ^ / dej] { - ^^(6) - (6) + Ti^(6)[ji^^ + Ji^^]}[/L(6) - fn{e)] 

r roo 

iLi = -r / E ^ / dc.ip„(u;i)r™(e,c.i){7V„(aa;i)/i(6)[l - 

^ a '.=±1 ^" (E3) 
+ 7V„(-a^i)/fl(e,„)[l-/L(e)]}, 

where Caa = e + crwi . The elastic transraission is given by 

ro(e) = Tr[G\e)Tn{e)G\e)TL{e)] (E4) 

and the inelastic prcfactor by 

T™(e,c.i) = Tr[G'(e,„)r«(e,„)G"(e,„)A"G"(e)ri(e)G'(e)A"], (E5) 
while the factors in the elastic correction are 

T-(6,c.i) = 2ReTr[G''(e)r;j(e)G"(e)ri(e)G''(e)A"G''(e,„)A"] 

T-^(e,c.i) - ImTr[G''(e)rfl(e)G"(e)ri(e)G''(e)A"G''(e,„)ri(e,„)G"(e,„)A"] (E6) 
T-«(e,c.i)=ImTr[G'\e)r^(e)G°(e)ri(e)G''(e)A"G''(e,„)r^(e,„)G"(e,„)A"]. 

The most complicated part are the integrals 
J!^{e) = j '^2Ke[Dl{u,)]ReTv[&'{e)TR{e)G\t)TL{e)&\e)\'^&\e~^^^^ 

(E7) 

and the coefficients coming from the last term in Eq. (|C7p for HJ' 

Ti\e) = 2ReTr[G\t)TR{e)G\e)TL{e)G\e)X'] 
Ji'''=Dl{Q) j ^Tr[G''(a;i)ri(^i)G"(a;i)A"]/iK) 



(E8) 



j'J'^^DliO) j ^Tr[G'^(a>i)ri^K)G°K)A"]/«(u;i). 



r 



In the equations above, the phonon density of states evaluated as principal part integrals. 
Pa{e) includes all possible broadening effects and shifts in the so-called wide band hmit, discussed in the text, 

of the bare vibrational frequencies. Without these effects the e and ui dependences of the coefficients Tq, T™, T^" 

Paie) = 6{e - fuva) - S{e + huja), and Eqs. ||E7|) must be and T^^^-^, may be dropped and they may be simply 
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evaluated at the Fermi energy. Then the integrals over 
products of Fermi functions in the corresponding current 
terms can be done analytically. One also finds that all 
terms of 51 ei which do not involve Eqs. (|E7|I or ljE8|l yield 
conductance contributions 5Gei (e) = dSIei / which are 
symmetric in the bias: 5Gei{~V) — 5Gei{V). We call the 
sum of these terms Sl^f"^ , and I'^v"^ = /o^ + Sl^f™^ + ■ 
The same may be done to the trace expressions in Eqs. 
(|E7|) . and one finds that the corresponding part in the 
current (J/gj yields the asymmetric current 7°'*^ [3^ . 



The contribution of Eqs. IjESp to the current is typically 
very small in the small-voltage limit which we are con- 
sidering. Furthermore, since they do not introduce any 
relation between the voltage and the vibrational frequen- 
cies, they cannot give a contribution to the conductance 
steps. Therefore, we drop them for the sake of simplicity. 
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